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High-order ADI schemes for 
convection-diffusion equations with mixed 
derivative terms 


B. Diiring, M. Fournie and A. Rigal 


Abstract We present new high-order Alternating Direction Implicit (ADI) schemes 
for the numerical solution of initial-boundary value problems for convection-diffusion 
equations with cross derivative terms. Our approach is based on the unconditionally 
stable ADI scheme proposed by Hundsdorfer ifT^ . Different numerical discretiza¬ 
tions which lead to schemes which are fourth-order accurate in space and second- 
order accurate in time are discussed. 


1.1 Introduction 

We consider the multi-dimensional convection-diffusion equation 


Uf = div(DVM) -fc • Vm 


( 1 . 1 ) 


on a rectangular domain £2 C supplemented with initial and boundary condi¬ 
tions. In (II.11 1. 
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are a given nonzero convection vector and a given, fully populated (non-diagonal), 
and positive definite diffusion matrix, respectively. Thus, both mixed derivative and 
convection terms are present in ( II.lb . 

After rearranging, problem (II. lb may be formulated as 


du(x,y,t) 

dt 


= {di2 + d2i) 


d^u du d^u du d^u. 

^ + {c,-+dn^) + {c2^+d22^). 


( 1 . 2 ) 


::Fo(«) 


= :Fi(«) 


=-Fi(u) 


This type of convection-diffusion equations with mixed derivatives arise fre¬ 
quently in many applications, e.g. in financial mathematics for option pricing in 
stochastic volatility models or in numerical mathematics when coordinate transfor¬ 
mations are applied. Such transformations are particularly useful to allow working 
on simple (rectangular) domains or on uniform grids (to have better accuracy). Thus, 
this approach allows to consider complex domains or to define non-uniform meshes 
to take into account the stiffness behavior of the solution in some part of the domain. 

In the mathematical literature, there exist a number of numerical approaches to 
approximate solutions to (fO . e.g. finite difference schemes, spectral methods, fi¬ 
nite volume and finite element methods. Here, we consider (II.lb on a rectangular 
domain Q. C In this situation a finite difference approach seems most straight¬ 
forward. 

The Alternating Direction Implicit (ADI) method introduced by Peaceman and 
Rachford m, Douglas ID |3, Fairweather and Mitchell El is a very powerful 
method that is especially useful for solving parabolic equations on rectangular do¬ 
mains. Beam and Warming El, however, have shown that no simple ADI scheme 
involving only discrete solutions at time levels n and /i + 1 can be second-order ac¬ 
curate in time in the presence of mixed derivatives (Fq 7 ^ 0 in (11.2b ). To overcome 
this limitation and construct an unconditionally stable ADI scheme of second order 
in time, a number of results have been given by Hundsdorfer 112 [HI and more re¬ 
cently by in ’t Hout and Welfert ifTOl . These schemes are second-order accurate in 
time and space. 

High-Order Compact (HOC) schemes (see, e.g. O [141) employ a nine-point 
computational stencil using the eight neighbouring points of the reference grid point 
only and show good numerical properties. Several papers consider the application 
of HOC schemes (fourth order accurate in space) for two-dimensional convection- 
diffusion problems with mixed derivatives EE but without ADI splitting. More¬ 
over, the HOC approach introduces a high algebraic complexity in the derivation of 
the scheme. 

We are interested in obtaining efficient, high-order ADI schemes, i.e. schemes 
which have a consistency order equal to two in time and to four in space, which are 
unconditionally stable and robust (no oscillations). We combine the second-order 
ADI splitting scheme presented in |[T2 [TOl with different high-order schemes to 
approximate Tb,Fi,F 2 in (11.2b . We note that some results on coupling HOC with 
ADI have been presented in |[T3l . however, without mixed derivative terms present 
in the equation. 
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Up to the knowledge of the authors there are currently no results for ADI-HOC 
in the presence of mixed derivative terms. In this preparatory work we validate the 
coupling of ADI and HOC by numerical experiments. 


1.2 Splitting in time 

In time, we consider the following splitting scheme presented in ifT^fTOl . We con¬ 
sider (ll.21 l. and we look for a (semi-discrete) approximation (/" ~ «(?„) with f„ = nA, 
for a time step A,. The scheme used corresponds to 

yi =yO-t-0A;(Ti(y')-T’i(f/"-i)), 
y2 = yi + 0A,(T2(y2) 

< y'> = yO-f (jA,(T’(y2)-T(i/"-')), (1.3) 

yi =y‘'4-0A,(Ti(y')-T’i(y2)), 
y2 = yi + 0A,(T2(y2)-T2(y2)), 

^t/" = y2, 

with constant parameters 0 and (7, and F = Fq + F\+ F 2 . To ensure second-order 
consistency in time we choose a — Xjl. The parameter 0 is arbitrary and typically 
fixed to 0 = 1/2. The choice of 0 is discussed in ifT^ . Larger 0 gives stronger 
damping of implicit terms and lower values return better accuracy (some numerical 
results for 0 = 1/2-1- v^/6 are given in section fU^ . 

We note that Fq is treated explicitly, whereas Fi,F 2 (unidirectional contributions 
in F) are treated implicitly. In the following section, we discuss different high-order 
(fourth order) strategies for the discretization in space. 


1.3 High-order approximation in space 

For the discretization in space, we replace the rectangular domain Q = x 

[L 2 ,R 2 ] C with Ri > L\, R 2 > L 2 by a uniform grid Z = {xi S : x,- = 

Li + (i - 1)A^, i = 1,. .. ,A} X {yj G [L 2 ,R 2 ] : yj = T 2 + U - i)^y, j = 
consisting of N xM grid points, with space steps Ax= {Ri — Li)/{N — 1) and Ay = 
{R 2 —L 2 ) / (M— 1). Let Hi j denote the approximate solution in {xi,yj) at some fixed 
time (we omit the superscript n to simplify the notation). 

We present different fourth-order schemes to approximateFb,Fi,F 2 in (IL3b . The 
first one uses five nodes in each direction and the second one is compact. Both 
schemes are considered with boundary conditions of either periodic or Dirichlet 
type. 
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1.3.1 Fourth-order scheme using five nodes 

We denote by 5x0, dx+ and the standard central, forward and backward finite 
difference operators, respectively. The second-order central difference operator is 
denoted by 5^, 




= 5.1 5, 




1 J '^^ij 1 J 

• 


The difference operators in the y-direction, 5vo, 5,,+, dy- and dy, are defined analo¬ 
gously. Then it is possible to define fourth-order approximations based on. 


{Ux)i.j 

{^xx)ij 

{^yy)i,j 


1 ^5,, ) 5yOM, j' = 

1-^52 ]5V-, = 


12 

c2 


^h'+2,7 “t” j' 8 m /—1 j' -f Ui—2J 

124^ ’ 

^i,2+2 “t” Suj j—i -f lti j—2 

i2a;^ ’ 

-Ui+2J + 16m,+ i - 30m, j -F 16m,_i j - m,_2,; 

1242 ’ 

-Ui,j+2 + 16m, j+i - 30 m, j -F 16 m,j_i - m, j_2 


1 


42 , 


1242 


1-^52 5 ^ 1 _ ^52 5 ^^,^ 


1444^4j, 


[64 (m,-^ 1 ,2+ 1 w,'— 1 ,2+ 1 “F M,_ \ j-_ \ M,-^ 1 1 ) 

-F8( —M,'+2,2+1 ~ “i+l.2+2 + «,-1.2+l + “!-2,2+1 

— “ i '-2,2-1 ~ “i-1,2-2 + “ i '+1,2-2 + “i+2,2-l) 
+ (“i+2.2+2 — “(-2.2+2 + «(-2.2-2 — «(+2.2-2)] ■ 


(1.4) 

For each differential operators appearing in Fq, F\ and F 2 , we use these five-points 
fourth-order difference formulae. 

Combining this spatial discretization with the time splitting ( ll.31 l. we obtain a 
high-order, five-points ADI scheme denoted H05. Its order of consistency is two in 
time and four in space. 


1.3.2 Fourth-order compact scheme 


We start by deriving a fourth-order HOC scheme for 


F\ (m) = d\ 


d^u 

dx^ 


-ci^ =g, 

ax 


(1.5) 


with some arbitrary right hand side g. We employ central difference operators to 
approximate the derivatives in (11.5b using 
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^{xi,yj) = dxauj 

Ar d^u, , /f. 

(1.6) 

d^u 2 

■^{xi,yj) = 8xUij 

A^ d^u, , ,/ts 

(1.7) 


By differentiating (11.5b . we can compute the following auxiliary relations for the 
derivatives appearing in ( 11.6b . ( 11.7b (in the following, for the sake of brevity we 
omit the argument {xi,yj) of the continuous functions) 

d^u 1 dg Cl d^u 

dx^ dll dx dll dx^ ’ 

1 d^g Cl d^u 1 d^g Cl f \ dg ci d^u 

dx^ dll dx^ dll dx^ dn dx^ dn yt/n (9x dn dx?- 

Hence, using (11.8b and ( 11.9b in (11.6b and (11.7b . respectively, equation (11.5b can be 
approximated by 



dii5^ Uij + ciSxoUij — gij + 


12 


( 

\r/ii dx 


d^g cj d^u\ 

dyp- dll dx?- J 


^(4^). (1.10) 


We note that all derivatives on the right hand side of (II.10b can be approximated 
on a compact stencil using second-order central difference operators. This yields a 
high-order compact scheme of fourth order for ( 11.5b which is given by 


diiSxUi j ClUi j + -- 1 - 

12 dll 




/ ~ SiJ + 


Cl 

12 V'^ii 


^xOgi,, 


c 8i,j 


( 1 . 11 ) 


In a similar fashion we can discretize the operator T 2 (m ) = g by a high-order compact 
scheme of fourth order given by 


dildyU^ 


cj 


■ CldyilUi^j 


12d22 


y 


= g: 


I,] 






( 1 . 12 ) 


Defining vectors U = (ui.i, ■ ■ -.um^m) and G = (gi.i,... ,gNM), we can state these 
schemes (II.11b and (11.12b in matrix form A^U = B^G (for Fi iu) = g) and AyU = 
ByG (for F 2 {u) = g), respectively. We apply these HOC schemes to find the uni¬ 
directional contributions T*, T*, and Y^, in ( fT3b . respectively. For example, to 
compute 

Y^ =Y^ + ^[Fi{Y^)-Fi{U''-^)) 

in the second step of (11.3b (which is equivalent to 721 ( 7 ^ — = — ^(7** — 7*)), 

we use Aj:(7* — = Bxi—^(Y^ — 7^)) that can be rewrite into 



= 47 " 
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Note that the matrix {B^ — {At/2)Ax) appears twice in ( 11.31 ). in steps two and five. 
Similarly, [By — [At/2)Ay) appears in steps three and six of ( 11.31 ). Hence, using LU- 
factorisation, only two matrix inversions are necessary in each time step of ( 11.3b . 
Moreover, for the case of constant coefficients, these matrices can be LU-factorized 
before iterating in time to obtain an even more efficient algorithm. 

To compute and in steps one and four of (11.3b which require evaluation of 
Fq (mixed term) we use an explicit approximation using the five-points fourth-order 
formulae ( 11.4b . 

Combining this spatial discretization with the time splitting (11.3b . we obtain a 
high-order compact ADI scheme denoted HOC. Its order of consistency is two in 
time and four in space. 


1.4 Numerical experiments 

We present numerical experiments on a square domain Q = [0,1] x [0,1] for two 
types of boundary conditions, periodic and Dirichlet type. The initial condition is 
given at time 7b = 0 and the solution is computed at the final time Tf =0.1 with 
different meshes A^ = Ay = h and different time steps Af. In our numerical tests we 
focus on the errors with respect to time and to space. 

In the first part, we consider the periodic boundary value problem considered in 
ifTOl . We implement the scheme detailed in ifTol based on second-order finite differ¬ 
ence approximations (referred to as CDS below) and compare its behaviour to our 
new schemes H05 (section [OT]) and HOC (section [T32|(. In the second part, we 
consider Dirichlet boundary conditions and restrict our study to the more interesting 
HOC scheme. In that part, we extend the splitting scheme to a convection-diffusion 
equation with source term. 


1.4.1 Periodic boundary conditions 

The problem given in Qo) is formulated on the domain Q. = [0,1] x [0,1]. The 
solution u satisfies (fTTT ) where 

c = D = 0.025 (]2), 

with periodic boundary conditions and initial condition u[x^y, Tq) = 

We employ the splitting (11.3b with <7=1/2 and 0 = 1/2. 

We first present a numerical study to compute the order of convergence in time of 
the schemes CDS, H05 and HOC. Asymptotically, we expect the error e to converge 

e = CA™ 


as 
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at some rate m with C representing a constant. This implies 
log(e) =log(C)+mlog(4,). 

Hence, the double-logarithmic plot e against A, should be asymptotic to a straight 
line with slope m that corresponds to the order of convergence in time of the scheme. 
We denote by £2 and £00 the errors in the / 2 -norm and /„o-norm, respectively. We 
refer to Table 1 1.1] for the order of convergence in time computed for different fixed 
mesh widths h € {0.1,0.0.025,0.00625} and time steps 4, € [7y/30,7y/90]. The 
solution computed for 4, = 7}/100 is considered as reference solution to compute 
the errors. The global errors for the splitting behave like C{Ai)'^. We also observe 
that the constant C only depends weakly on the spatial mesh widths h. 


Table 1.1 Numerical convergence rates in time for 9 = } 


Scheme / 2 -error convergence rate /„o-error convergence rate 

h = 0.\ /7 = 0.025 /j = 0.00625 h = 0.l /; = 0.025 /? = 0.00625 


CDS 

2.2002 

2.1975 

2.1969 

2.1973 

2.1958 

2.1956 

H05 

2.1999 

2.1973 

2.1969 

2.1992 

2.1953 

2.1955 

HOC 

2.2002 

2.1973 

2.1969 

2.2007 

2.1953 

2.1955 


In the following, we study the spatial convergence. The double-logarithmic plots 
£2 and £00 against h give the rates of convergence. Contrary to the time convergence, 
the order now depends on the parabolic mesh ratio = Atj so the numerical 
tests are performed for a set of different constant values of /i. For simulations, ji is 
fixed at constant values jx G {0.4,0.2,0.1,0.005} while Ax = Ay = h^Q {At is then 
given by At = lih^). The results for the / 2 -error are given in Table [T2l and for the 
/ 00 -error in Table fOl The solution computed for h = 0.00625 is used as reference 
solution to compute the errors. 


Table 1.2 Numerical convergence rates in space of / 2 -error for fixed fi as Ax, A, —> 0 and 9 = } 


Scheme 

p = 0.4 

p = 0.2 

p = 0.1 

p =0.05 

CDS 

1.7828 

1.7909 

1.7821 

1.7845 

H05 

2.2291 

2.5188 

2.8153 

3.0672 

HOC 

2.2685 

2.5191 

2.8152 

3.0671 


Remark: The choice of the parameter 0 is discussed in ifT^ . However, for the 
convergence rates, 0 seems to have little influence. For example, for the scheme 
H05 with 0 = 1/2 -b v^/6 we obtain very similar results as shown in Table 1 1.4] 
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Table 1.3 Numerical convergence rates in space of /^-error for fixed ji as Ax, A, —> 0 and 6 = 5 


Scheme 

p = 0.4 

p = 0.2 

p = 0.1 

p = 0.05 

CDS 

1.7170 

1.7125 

1.7040 

1.7038 

H05 

2.2931 

2.6166 

2.9182 

3.1584 

HOC 

2.3175 

2.6176 

2.9184 

3.1584 

Table 1.4 

Numerical convergence rates in space for H05 for fixed fl as Ax, A, —> 0 and 9 


p = 0.4 

p = 0.2 

p = 0.1 

p = 0.05 

I 2 rate 

2.2310 

2.5186 

2.8152 

3.0671 

/oo rate 

2.2938 

2.6164 

2.9181 

3.1584 


1.4.2 Dirichlet Boundary conditions 

In this section we only consider the HOC scheme which presents more interesting 
properties than the other schemes. Indeed, compared to CDS, its accuracy is larger 
and compared to H05, no specific treatment at the boundaries is required for the uni¬ 
directional terms Fi, F2, the compact scheme is optimal in this respect. A particular 
treatment is necessary when ghost points appear in the explicit approximation of 
the mixed term Fq. To preserve the global performance, the accuracy of the approx¬ 
imation near the boundary conditions has to be sufficiently high. We have used a 
sixth-order approximation in one direction (although lower order may also be used 
191). For example, for uq j on the boundary, at a ghost point m_i j we impose 

“-1.; = 5 mo.; ~ IOmi,/ + 10m2j' — 5 uT,.j + M4,;. 

For the numerical tests, we consider the problem 

lit = div(DVM) +c-'Vu + S 
on the domain Q = [0,1] x [0,1] where 

" = -(]). 0 = »» 25 (] 4 )- 

and the source term S is determined in such a way that the solution is equal to 
u{x,y,t) = — sin(a:jc) sin(7ry). The Dirichlet boundary condition and initial con¬ 
dition are deduced from the solution. To incorporate the source term S in the splitting 
(fO . F needs to be replaced by F + 5'. More specifically, F[U" is replaced by 
F -\-S{t''^^) andF(T^) by F{Y'^) + S[F). We perform the same numerical ex¬ 

periments as in the previous section. The final time is fixed to 7y = 0.1 and the errors 
are computed with respect to a reference solution computed on a fine grid in space 
(Ax = Ay = 0.00625). Different meshes in space are considered for Ax = Ay = h and 
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h £ {0.1,0.05,0.025,0.0125}. For pi = 0.4 the double-logarithmic plots £2 and & 
against h are given in Figure [TTI 




Fig. 1.1 Numerical convergence rate in space for HOC (9 = }) and p = 0.4. 


The results of several numerical tests are reported in Table ll.5l for fixed parabolic 
mesh ratio ji = AtfA^ while Ax, At —0. In all situations, the new HOC scheme 
shows a good performance with fourth-order convergence rates in space, indepen¬ 
dent of the parabolic mesh ratio fi. 


Table 1.5 Numerical convergence rates of / 2 -en'or and /^o-en'or for HOC (9 = }) for different 
constant values of p (Dirichlet boundary conditions). 



p = 0.4 

p = 0.2 

p = 0.1 

p =0.05 

I 2 rate 

4.0971 

4.1875 

4.2129 

4.2196 

U rate 

4.1530 

4.2372 

4.2717 

4.2806 


1.5 Conclusion 

We have presented new high-order Alternating Direction Implicit (ADI) schemes for 
the numerical solution of initial-boundary value problems for convection-diffusion 
equations with mixed derivative terms. Using the unconditionally stable ADI scheme 
from m we have proposed different spatial discretizations which lead to schemes 
which are fourth-order accurate in space and second-order accurate in time. 

We have performed a numerical convergence analysis with periodic and Dirich¬ 
let boundary conditions where high-order convergence is observed. In some cases, 
the order depends on the parabolic mesh ratio. More detailed discussions of these 
schemes including this dependence and a stability analysis will be presented in a 
forthcoming paper. 
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